Critical jamming of frictional grains in the generalized isostaticity picture 
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While frictionless spheres at jamming are isostatic, frictional spheres at jamming are not. As a 
result, frictional spheres near jamming do not necessarily exhibit an excess of soft modes. However, a 
generalized form of isostaticity can be introduced if fully mobilized contacts at the Coulomb friction 
threshold are considered as slipping contacts. We show here that, in this framework, the vibrational 
density of states (DOS) of frictional discs exhibits a plateau when the generalized isostaticity line is 
approached. The crossover frequency uj* scales linearly with the distance from this line. Moreover, 
we show that the frictionless limit fi — » 0, which appears singular when fully mobilized contacts are 
treated elastically, becomes smooth when fully mobilized contacts are allowed to slip. 

PACS numbers: 45.70.-n, 46.55.+d, 63.50.-x, 81.05.Rm 



The jamming transition in disordered media has been 
the focus of intensive research efforts in recent years. For 
frictionless spheres, the transition occurs at the isostatic 
point where the contact number z reaches z? = 2d 
This point, known as point J, acts as a critical point: 
close to it the spectrum of vibrational modes shows an 
excess density of states (DOS) at low frequencies, with 
a finite number of zero-energy modes at the jamming 
transition and a cross-over frequency u>* ~ (z—z? so ) 

Friction has to be introduced to make connection with 
the granular systems studied in experiments 0], in en- 
gineering [f| and in the geosciences 0]. For frictional 
spheres, the contact number at the transition lies in the 
range z? — d+1 < z < 2d, depending on the preparation 
method and the friction coefficient \x of the material. The 
frictional isostatic value z^ so derives from a counting ar- 
gument assuming arbitrary tangential friction forces, and 
is reached in practice at jamming in the limit fi — > oo, 
for gently prepared packings 0, d, • 

However, the magnitude of the tangential forces is lim- 
ited by the Coulomb criterion |/ t | < (if n - Let m — 
\ft\l 'nf n be the mobilization of a contact, such that con- 
tacts at the Coulomb threshold, so-called fully mobilized 
contacts, have m = 1. Crucially, there can be a finite frac- 
tion of contacts at the mobilization threshold, i.e. with 
fixed ratio \ft\/nfn 11 1- This affects the counting 
arguments, and suggests to consider a generalized iso- 
staticity condition [7|. 

We define n m £ [0, z/2] as the fraction of contacts per 
particle at the Coulomb threshold. For frictional par- 
ticles, the tangential forces introduce d — 1 additional 
force components at each contact. Then for a packing 
to be stable, the Nd(d + l)/2 rotational and transla- 
tional degrees of freedom need to be constrained by the 
Nzd/2 — Nn m independent force components. This lead 
Shundyak et al. 0] to propose the generalized isostaticity 
criterion 

z>(d+l) + ^=zZ- (1) 
Simulations have shown that as frictional packings are 




FIG. 1: Position of our packing configurations in the phase 
space defined by z and n m in our 2d simulations. The gen- 
eralized isostaticity line n m = d(z — z™ )/2 — z — 3 is shown 
in black. Stable packings can not exist left of this line. Each 
marker denotes one of 30 configurations at a given (u, p) with 
the color-marker combinations indicated by the legend. For a 
given fi (a given symbol), configurations with larger p are to 
the right, and n m drops to zero for the largest pressures. 



prepared gently, the fraction of fully mobilized contacts 
at jamming indeed is such that z approaches the general- 
ized isostaticity line z = z™ in the z-n m plane (Fig. [!}, 
in accord with an earlier suggestion by Bouchaud [l(j |. 

Fully mobilized contacts can not resist tangential per- 
turbations, and simulations evidence that a substantial 
part of the movement of frictional piles is due to contacts 
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failing at the mobilization threshold [Bl 11. 12. 
exact values of for < [i < oo depend on the prepa- 
ration method [{j, [l2| , but in the limit of very "gentle" 
equilibration so as to approach the jamming threshold, 
they are a smooth decreasing function of \i @, H, . For 
2d systems, one finds that n m —>■ 1 in the frictionless limit 
since z « z® so = 4, while n m — > in the limit [i — > oo, 
as z then approaches the frictional isostaticity criterion 
z = d+1. 

In this paper, we show that the similarities between 
frictionless and frictional static sphere packings, which 
are brought to the foreground by this concept of general- 
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ized isostaticity, also extend to the dynamic properties. 
We do so by calculating the density of states (DOS) of 
frictional sphere packings, while taking into account that 
contacts at the mobilization threshold m = 1 slip with un- 
changed tangential forces during small amplitude vibra- 
tions. We recover the low- frequency plateau in the DOS 
which is characteristic for packings of frictionless parti- 
cles near jamming. We find that the rescaled crossover 
frequency uj* = uj* /p 1 ^ 6 scales linearly with the distance 
from the generalized isostaticity line, as ui* ~ (z — z™ ), 
see Fig. [5^. An analysis of the eigenmodes shows that 
the dominant low-energy displacements consist of parti- 
cles rolling without sliding for contacts with m ^ 1, and 
tangential sliding at contacts with m = 1. 

This scenario nicely generalizes the findings for fric- 
tionless spheres to the frictional case. Our findings stress 
the crucial role of the response of fully mobilized con- 
tacts. When instead, these are all taken to be elastic, 
the inclusion of friction becomes a singular perturbation, 
and the excess low-frequency modes seen in frictionless 
packings are suppressed 
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Simulation method — We perform a numerical analy- 
sis of two-dimensional packings of frictional particles in- 
teracting with Hertz-Mindlin forces (during the prepara- 
tion of the packings, energy is dissipated through viscous 
damping — see [7J, LLJ]). The Hertz-Mindlin interaction 
is comprised of a normal interaction that scales as 5 3 ^ 2 , 
where 6 is the overlap of the particles, and an incremental 
tangential interaction which is limited by the Coulomb 
criterion f t < pf n (15j |. The packings are made at fixed 
friction coefficient fi and pressure p. 

Motivated by ([T]), we map the configurations in the 
phase space of Fig. [T] spanned by z and n m . The gener- 
alized isostaticity criterion z — z™ then defines a line of 
marginal packings with end points [z = 3, n m — 0) and 
(z = 4,n m = 1), so that stable packings lie to the right 
of it. As in [7], we observe that for sufficiently gentle 
preparation methods, packings approach the generalized 
isostaticity line in the limit p — > 0. Hence we identify 
this line with the (zero pressure) jamming transition for 
slowly equilibrated frictional discs. For larger pressures, 
packings attain higher contact numbers, and we find a 
range of pressures where the distribution of n m is bimodal 
(see e.g. the configurations for (fi = 0.001, p = 2.4- 10~ 4 ), 
due to the interplay of elastic and viscous forces during 
the preparation of the packings Q . We focus here on the 
behavior for p — * at fixed viscous damping. 

Treatment of small amplitude vibrations — The vi- 
brational density of states is obtained by linearizing the 
equations of motion around each stable state. Rattlers, 
which give rise to trival zero-frequency modes, are left 
out from the analysis from the start. 

There are several subtleties associated with frictional 
response. First of all, the nature of the ideal Coulomb 
condition |/ t | < \if n implies a discontinous response for 
sliding displacements at fully mobilized contacts. For 




FIG. 2: (a) Normalized crossover frequency uj* where the 
plateau in the DOS is reached as a function of the distance 
z — (3 + n m ) from the line of generalized isostaticity. The 
relation is linear, in spite of the change of the /i-n m -relation 
at larger pressures (see Fig. [1} . (b) Contact geometry illus- 
trating the various displacements and rotations as well as the 
effective tangential displacement St of the contact point. For 
simplicity, we have assumed that particle i is stationary. 



displacements which lead to an increase in the tangential 
force ft (taking /„ fixed), the Coulomb condition then 
implies that the contact slips with ft unchanged — this 
translates into a contact stiffness k t = for this dis- 
placement. For displacements in the opposite direction, 
however, f t will decrease so that kt is effectively nonzero! 

This separation into two types of displacements at 
fully mobilized contacts is only meaningful for static re- 
sponse studies. Under vibrations, modes effectively cou- 
ple through the Coulomb condition; moreover fully mo- 
bilized contacts dissipate energy during the slipping half 
of the phase, and a decomposition into purely nondissi- 
pative eigenmodes is not possible. In order to avoid these 
complications — which are artifacts of the singular na- 
ture of the ideal Coulomb condition — we here simply 
put the tangential stiffness fc t = for all fully mobilized 
contacts: fully mobilized contacts are treated as slipping 
contacts in both directions. This allows us to study the 
DOS of these modes. 

Secondly, a detailed analysis of the dynamics shows 
that there always is, even at contacts which are not fully 
mobilized, a nonpotcntial term in the dynamical equa- 
tions due to the tangential friction Q . The origin of this 
term is the change in the direction of the tangential fric- 
tion force when a contact point rolls over a particle. This 
effect is of the order of ft itself, which for our Mindlin 
forces scales as ~ 8 3 ^ 2 . Hence for packings close to jam- 
ming, the contribution of these terms can be neglected in 
comparison with the normal and tangential springs which 
scale as 5 1 / 2 ; for this reason we will disregard these terms 
in our analysis below (similar to what is often done with 
the pre-stress terms in frictionless packings). 

With these approximations, the equations of motion 
are conservative to first order, and their structure re- 
sembles what one finds for the vibrational properties of 
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frictionless particles. After we expand the equation of 
motion around the equilibrium we obtain equations of 
motion of the form Sr a = —D a p5rp + 0{8r 2 ), where 
the dynamical matrix D a p can, with the present simpli- 
fications and after a rescaling of the coordinates by the 
square root of the masses/moments of inertia, be writ- 
ten in terms of the derivatives of an effective potential, 
D aP = {m a m )~ 1/2 ds f X js[. For our packings with 



Hcrtz-Mindlin forces the effective potential is given by 



V =^2 K(5r.nf ~ ^(Sr.t) 2 + k t 6t 2 



(2) 



(ij) 



where St = (Sr.t) — (RiSai+RjSctj) is the tangential dis- 
placement of the contact point of the two particles, as 
illustrated in Fig. [2fb) . Here k n and kt are the normal 
and tangential stiffness, respectively, which derive from 
the Hertz-Mindlin interaction law and which both scale 
as In this formulation, the vibrational 

spectrum is a histogram of the eigenvalues lo\ of D a p as 
a function of the associated frequencies lo^. 

Density of States — We have studied the vibrational 
density of states (DOS) for a wide range of friction coeffi- 
cients (p, G [1CT 3 , IP 3 ]) and pressures (p <E [10~ 6 , 10~ 3 ]). 
As in previous work [14J, we plot the density of states 
as a function of the rescaled frequencies Q — to/p 1 / 6 , ap- 
propriate for Hertz-Mindlin forces which exhibit a trivial 
softening with frequencies scaling as \fk ~ p 1 / 6 . 

Figs. (3^-b show the DOS for the smallest pressure 
(p = 1.41 • 10 -6 ), i.e, close to the generalized isostaticity 
line, for the full range of /x. When the fully mobilized 
contacts are allowed to slip, the density of states clearly 
shows an excess number of low-energy modes (Fig. 
Otherwise, when all contacts are chosen to have a finite 
kt, the plateau in the DOS disappears for small friction 
coefficients, where n m is large (Fig. [3)3 and [141]). 

In Fig. [3b, we show the evolution of the density of 
states, again allowing fully mobilized contacts to slip, 
with increasing pressure for \i — 0.3. The low-frequency 
part of the DOS shows linear D(Cu) ~ Co Debye-like be- 
havior up to a normalized frequency Co* which increases 
with pressure. In the frictionless case, Co* marks the fre- 
quency scale below which the packing can be treated like 
an elastic solid and above which the DOS shows a plateau 
For our frictional packings, an extension of this ar- 
gument then predicts the scaling Co* ^ z — z^ Q with in 2d 
the generalized isostatic contact number z ; ™ = 3 + n m 
— we will verify this prediction below. 

We extract Co* in the following way: Since the DOS at 
different pressures do not have similar functional forms, 
we cannot rescale the integrated DOS as was done in [14 1. 
Instead, we smooth the integrated DOS to obtain an in- 
terpolated DOS. Then Co* is determined by the point at 
which D(Co) reaches a value of 0.2, normalized to the 
height of the plateau at Q = 1, to avoid nonlinearities in 
the approach to the plateau. Fig. [5b, shows the Co* we 




FIG. 3: Density of states (DOS), (a) DOS with fully mobilized 
contacts treated as slipping for the smallest pressure p — 1.41- 
10~ 6 (approaching the line of generalized isostaticity) for a 
range of p.. (b) DOS for the same packings as in (a), but with 
all contacts treated as non-slipping, as in [l4| . (c) Illustration 
of the low frequency plateau developing in the DOS as p is 
decreased, for an intermediate friction coefficient p = 0.3 and 
with m = 1-contacts slipping, (d) DOS for range of pressures 
and \x = 0.001, with the tangential stiffness constant set to 
(see text). 
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obtain for different \x and p as a function of z — z™ . The 
relation is linear to a good approximation, confirming our 
prediction. 

Even if the fully mobilized contacts are allowed to slip, 
the density of states for fj, — > is still noticeably different 
from the frictionless case with fi = (although both have 
an excess of low-density modes) . This is because the non- 
mobilized contacts still have a finite tangential stiffness 
kt comparable to the stiffness for compression of bonds 
k n , and hence a finite influence — this is even true when p 
approaches zero and the system approaches generalized 
isostaticity. Clearly, this non-smooth behavior has its 
root in the singular change of the dynamical matrix, due 
to the finite value of k t of many contacts 1J|. Fig. [3]i 
illustrates that when one takes the limit k t — * in the 
dynamical matrix, one recovers the frictionless DOS, with 
a (5-function of weight N due to trivial rotational modes 
at Q = [13. 

Nature of the low- energy displacements — We now in- 
vestigate the nature of the eigenmodes of the low-energy 
eigenvalues of the dynamical matrix. Eq. @ predicts 
the nature of the lowest-energy displacements in the limit 
Co — > 0, p — * 0. The prefactors of the tangential, normal 
and sliding displacements scale as <5 3 / 2 , S 1 ^ 2 and S 1 ^ 2 , re- 
spectively, for a contact with m/1. For a m = 1-contact, 
kt = 0. Therefore, in the limit discussed above, the only 
low-energy displacement allowed for m =^ 1 is a purely 
tangential motion in combination with rotations of the 
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FIG. 4: Sample displacement for a low-energy eigenvector at 
p = 1.41 • 10 -6 and fj, = 0.001. (a) Contact network in the 
initial state, with fully mobilized contacts in red. (b) Parti- 
cles have been displaced proportional to the mode amplitude. 
The lines now link the particle centers to the position of the 
former contact points after rotation. If St ^ 0, the contact 
line between particles is broken. 



particles in such a way that there is no tangential sliding. 
For m = 1 contacts, the movement still has to be tan- 
gential, but sliding at the contact is permitted. A simple 
illustration of the displacements in a low-energy, low p 
mode is shown in Fig. 2J the difference in the rotational 
response between slipping fully mobilized contacts and 
non-slipping non-mobilized contacts is clearly visible. 

The occurrence of floppy modes at generalized iso- 
staticity, and their spatial structure which locally ex- 
hibits rolling and sliding depending on the contacts mo- 
bilization, can be understood from a counting argument. 
The 37V degrees of freedom of the packing participating 
in a floppy deformation need to satisfy two groups of 
constraints: setting the normal motion at each contact 
to zero gives Nz/2 constraints, while setting the sliding 
motion at each m ^ 1-contact to zero gives Nz/2 — Nn m 
constraints. Hence the number of constraints of the mo- 
tion is exactly equal to the number of degrees of freedom 
only if z = z™ Q = 3 + n m — floppy modes arise at gener- 
alized isostaticity. Note that the variational argument of 
Wyart et al. Q that estimates Co* in terms of distorted 
floppy modes then generalizes in our case to to* ~ z — z™ a . 

In the limit fj, — ► 0, we observe that n m — > 1 (every 
other contact is fully mobilized, see also [7]), so that 
at that point on the generalized isostaticity line, the 
Nz/2 — Nn m = N contacts with 1 equal precisely 
the number of rotational degrees of freedoms. Hence 
there are precisely enough contacts to couple the rota- 
tions to the translations. Unlike in the case of the spher- 
ical limit of ellipsoids, where the rotational modes appear 
in the gap of translational modes [It], EH , here transla- 
tions and rotations remain strongly coupled. 

Conclusion and outlook — The jamming transition of 
frictional packing of spheres is more complex than the 
transition at point J for frictionless spheres. However, 
we show that along the generalized isostatic line in the 
space spanned by the contact number z and the number 
of fully mobilized contacts n m , the system shows crit- 



ical behavior similar to the frictionless case if slipping 
contacts at m = 1 are incorporated into the dynamical 
matrix. In this case, the DOS shows a plateau near the 
z™ -line, and at larger pressure, the crossover frequency 
scales as Co* ~ (z — z™ ). 

The relation between the friction coefficient and n m 
is poorly understood except for the limiting cases fi — > 
oo and fi — > 0. Our simulations indicate that this rela- 
tion depends significantly on the preparation method of 
the sample, where packings prepared with higher viscous 
damping have larger n m [9J. The z-n m phase diagram 
can be a tool to understand the behavior of more com- 
plex packings. For example, for a packing undergoing 
an avalanche, the number of fully mobilized contacts in- 
creases before the system rearranges, so that the system 
moves along a vertical path in the z-n m space which even- 
tually crosses the z-^-line [12]. 
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